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Abstract 


A new computational technique for the solution of the full 
potential equation is presented. The method consists of outer 
and inner iterations. The outer iterate is based on a Newton like 
algorithm, and a preconditioned Minimal Residual method is used 
to seek an approximate solution of the system of linear equations 
arising at each inner iterate. The present iterative scheme is 
iormulated so that the uncertainties and difficulties associated 
with many iterative techniques, namely the requirements of 
acceleration parameters and the treatment of additional boundary 
conditions for the intermediate variables, are eliminated. 
Numerical experiments based on the new method for transonic 
potential flows around NACA 0012 airfoil at different Mach 
numbers and different angles of attack are presented, and these 
results are compared with those obtained by the Approximate 
Factorization technique. Extention to three-dimensional flow 
calculations and application in finite element methods for fluid 
dynamics problems by the present method are also discussed. 

I. Introduction 

The ability to compute transonic flow fields around airfoils 
or wings is an important aid in the design of efficient modern 
transport aircrafts since they operate predominantly in transonic 
ranges. Considerable effort has been spent, in recent years, on 
the construction of fast and accurate numerical procedures for 
the solution of the full potential equation. To be useful as a 
design and analysis tool, the success of a computational 


procedure should not be problem dependent. For example, some 
numerical procedures yield rapidly converged solutions if optimal 
values of acceleration parameters are provided and if other 
special conditions are given. However, it should be pointed out 
that optimal values of these parameters are often unobtainable 
for practical calculations. 

The standard iterative procedure for transonic small 
perturbation and full potential calculations was based on the 
successive line over-relaxation (SLOR) method. Because of its 
s^low convergence rates for many practical problems, the method 
has been replaced by many new iterative procedures. One of the 
most successful numerical techniques is based on the Approximate 
Factorization (AF) scheme, and there are many variants of the AF 
method^"* including those based on ADI ^ type developed by 
Ballhaus and Steger, AF2^ type by Ballhaus, et al., AF3‘ type by 
Baker, and SIP® type by Sankar, et al. These computational 
procedures provide substantial improvement in rates of 
convergence compared to the SLOR method. However, they all 
require one or more iteration parameters in order to accelerate 
the convergence, and an intermediate variable is also introduced 
into the iterative process for a two-dimensional calculation. 
Consequently, the uncertainty as to what values should be used 
for the iteration parameters, and the uncertainty about how to 
select the boundary conditions for the intermediate variable, may 
affect the convergence rates as well as the stability of the 
iterative process. It is our aim here to present an efficient 
iterative procedure which yields a rapid rate of convergence like 
the AF scheme, while eliminating the difficulties associated with 
the AF scheme. The present method consists of outer and inner 
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iterations. The outer iterate is based on a Newton like 
iterative process in which the Jacobian matrix is not required, 
and a preconditioned Minimal Residual algorithm is applied only 
to seek an approximate solution of the system of linear equations 
arising at each inner iterate. This method can therefore be 
regarded as a Newton like - Minimal Residual algorithm or an 
Inexact Newton like (IN) iterative procedure. 

The idea of the IN iterative scheme was first proposed by 
the author in [7], Although our early paper indicated that the 
method can be used to compute transonic flow fields around 
airfoils, it was not competitive with the AF scheme implemented 
by Dougherty et al®. The computational results showed that more 
iterations were needed for a converged solution compared to the 
AF scheme, and the CPU time per iteration for the IN method was 
about three times that required for the AF scheme. More 
recently, the IN method has been modified to include a better 
preconditioning operator so that a substantial improvement has 
been achieved: the number of iterations is now about half of 
that required by the AF scheme, and the CPU time per iteration is 
about twice of that required by the AF scheme. The present paper 
is mainly concerned with the numerical solutions of a 
two-dimensional full potential equation, and particular attention 
is focused on the improved version of the IN iterative method. 
Comparisons of numerical results for lifting and non-lifting 
airfoil calculations between the IN and AF schemes are given and 
extension to three-dimensional calculations will be discussed. 
We describe the problem formulation for the transonic flow fields 
around an airfoil in section II, and the solution of nonlinear 
discrete potential equations by the IN method is presented in 
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section III, results of computational experiments are discussed 
in section IV, and finally, concluding remarks are given in 
section V. 


II. Problem Formulation 

For a two-dimensional problem in Cartesian coordinates, the 
governing partial differential equation for an inviscid 
isentropic fluid flow expressed in the conservation form is 


(H ) + (H ) = 0 

XX y y 

where 

y - 1 2 2 l/(y-l) 

/> = [ 1 - (<}) + <t> )] 

y + 1 X y 


Equation (1) is known as the full potential equation, where ^ is 
the velocity potential, />, the density of the fluid flow, and y, 
the ratio of specific heats. Equation (1) is a nonlinear equa- 
tion since /> is a function of and <f> . The numerical solu- 

X y 

tion of Equation (1) for transonic flow is more delicate and more 
interesting than those for purely subsonic or purely supersonic, 
because the governing equation changes its type from elliptic in 
subsoinc regions to hyperbolic type in supersonic regions, and 
the boundary between these regions is unknown. Moreover, the 
equation also admits discontinuous solutions, such as shocks 
which may exist in the flow fields. 

To handle a general flow problem with complex geometries it 
is advantageous to transform Equation (1) from the physical 
domain in the Cartesian coordinates into the computational domain 
in a rectangle*. The full potential equation written in the 
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computational coordinates ^ and n is given by 


PU fiV 

( — ) + ( — ) = 0 ( 2 ) 
J t J n 


where 


r -1 

P = [1 

r +1 

U = A 4> + A <t> 

1 « 2 T, 


J = i n ~ ^ n 

X y y X 


l/(r-l) 

(U<f> + V(}> ) ] 
t n 

V = A <}> + A <t> 

2 ^ 3 n 

2 2 
A = t ^ 

1 X y 

2 2 


A=^n~tn , A=n + n 
2 xx yy 3 xy 

Here U and V are the contravar iant velocity components along the 
t and n directions, J is the Jacobian of the coordinate 
transformation, A^, Aj and A 3 are metric quantities. 

One of the difficulties in the numerical solution of 
transonic flow calculations is that both compression and 
expansion shocks are admitted by Equation (1). The expansion 
shocks, however, are physically meaningless. Thus in order to 
eliminate the expansion shocks from the flow fields, an 
artificial viscosity term is introduced, via an upwind bias, into 
the full potential equation. In this paper, the method of 
artificial density®"’ is implemented, where the fluid density is 
modified in such a way that 


A (/> - mA At) 


(3) 


where 

n = max [0, 1 - 1/M* ] 

Here s (t) indicates that s is replaced by t. In the above 
expression is a switching function which is zero in subsonic 
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flow fields and non-zero in supersonic flow fields, M- is the 
local Mach number, and /> is the density gradient in the upwind 
diretion. An important advantage in using the artificial 
density method is that a central difference approximation can be 
employed to discretize the full potential equation in the entire 
flow fields regardless of whether it is in a subsonic or a 
supersonic region. 


III. Solution Procedure 

By the application of the finite difference method the 
solution of the full potential equation (2) is transformed to the 
solution of a large set of nonlinear equations 

L(<|)) = 0 (4) 
where <|> is a vector of velocity potential at the grid points, and 
L is the nonlinear full potential operator. 

Newton Like Algorithm 

Our iterative scheme for the solution of Equations (4) can 
be described as follows: 

Let (})° be an initial guess for the velocity potential 
vector, compute the residual vector r° = L{<})°), then for 
n=0,l,2,..., until ||r "||2 < c. 


n n 


solve 

M 6<t> 

= -r 


n 



n + 1 

n 

set 

4 > 

= 


n+1 


compute 

r 

= L(4> 


where n is an iteration number, 6 ^ 
M„ is a matrix operator which varies 
It should be noted that if M„ = L' 


(5) 


is the correction vector, and 
from iteration to iteration, 
the Jacobian matrix of 
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L((|>), then (5) is a Newton's iterative process for the solution 
of the nonlinear equations (4). Although Newton’s method yields 
a rapid convergence rate, the method requires the initial guess 
<}>® to be inside a domain of attraction, that is, one must have a 
good initial vector to ensure for convergence. Furthermore, even 
if the linearized full potential operator is a sparse matrix, the 
Jacobian matrix L' (<|)) will likely be a full matrix. For many 
practical problems in aerospace industry the order of the 
nonlinear equations is large, such as 5000 or more, consequently 
it would be very difficult and expensive to compute the Jacobian 
matrix for each iterate n. 

In order to implement the iterative scheme in (5) 

efficiently it would then be natural to consider another operator 
for M„ . Axelsson^® has shown that if M„ is a linear operator and 
in some sense makes | |L( <}>") - | | almost insensitive to (}>", 

then the iterative procedure (5) converges. In this paper we 
shall choose M to be an approximation to the full potetntial 
operator, and with this particular choice the iterative process 
defines a Newton like algorithm. The construction for M„ is 
given as follows: 

Consider at the nth iterate, the fluid density has been 
calculated from values of the velocity potential at the (n-l)th 
iterate, the result of the application of a central difference 
approximation to L then leads to a nine-point formula, where 

(L<f>), =C (}) +W <|> +E (f) 

i/j i/j i/j i,j i-l/j i,j i+l,j 

+ N(j> +S(f> +NW4> (6) 

i/j i/j+1 i/j i/j-1 i,j i-l,j+l 

+ NE <|> + SW <|) + SE (J> 

i,j i+l,j+l i,j i-l,j-l i,j i+l,j-l 
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Note that, the values at the NW, NE, SE and SW positions are 

I 

usually much smaller than those at N, W, C, E and S positions, 
since they are due to the skewness effect of the coordinate 
transformation. The operator M is now chosen by setting the 
values at NW, NE, SW and ' SE to zero. Hence M will be a 
five-point formula and this i|'mplies that the skewness effect has 
been ignored. We would like to remark that neglecting the 
non-zero values at NW, NE,i SW and SE positions for the full 
potential equation in non-conservation form will no longer 
g^uarantee that M will be a good approximation to the L operator. 
Thus M should retain the nine-point formula structure for this 
situation, and the results for the non-conservative full 
potential equation will be reported later. 

It should be noted that for the conservative full potential 
equation M will be identical to the linearized operator L 
provided an orthogonal transformation is used, and for other 
transformations M will only be an .approximation to L. Thus the 
operator M takes the following form for orthogonal 

transformations : 

-e />Ai ^ ^ PA3 

M6<}> = [a ( — ) a + a ( — ) a ]6<t> (7) 

i/j € J i + l/2,j f: n J i,j+l/2 n -i,j 

in purely subsonic flow calculations. For mixed subsonic- 
supersonic flow problems, the density p in M has been modified 
according to (3) so that an artificial viscosity term is 
introduced. However, for large supersonic regions (i.e. strong 
shock calculations) it is necessary to introduce an additional 
upwind directional bias in the supersonic flow fields to ensure 
for smooth convergence. This can be achieved by modifying the 
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operator M so that a ^ type term 

it 

will take the form: 


is explicitly 


included, and M 


M6(f) = [± fi$d + d ( ) d 

i/j i e. J i + l/2,j t 


( 8 ) 


^ 3 -»■ 

+ d ( ) d ] 6 ^ 

n J i,j + l/2 Ti i,j 

for transonic flow calculations, where ^ is a switching function 
which is zero in subsonic regions and is non-zero in supersonic 
regions, ^ is a constant which controls the amount of the 

€t 

term that is introduced. 

^ ' 

It should be noted that Equations (7)-(8) are valid only for 
orthogonal transformations. In this paper a non-orthogonal grid 
transformation is used and the operator M will only be an 
approximation to L. Consequently, instead of Equation (8), we 
have: 


M6<j> = + L + EJ 6(|) (9) 

i, j ^ i, j 

where E is the error matrix due to ignoring the skewness effect 
in the grid transformation. Now consider E = E^ + al , then the 
Newton like iterative scheme in (5) becomes; 


n+1 n n 

i±H$d 4L + oI+E)(<|> “<t>) = ~L<}) (10) 

t 1 

Since l.|Ei|| << ||L|| for the full potential equation in 
conservation form, it is not hard to observe tht Equation (10) in 
fact simulates a time-dependent problem; 


n+1 

c<f) ± + L(}) =0 

t tt 

We would also like to remark that our iterative 


( 11 ) 


scheme is 
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fully implicit, and moreover, the boundary conditions for the M 
operator are the same as that imposed for the full potential 
operator. Although the iterative scheme given in (5) appears to 
be similar to that based on the AF technique, there is an 
important difference in choosing the operator for M. In the AF 
scheme M is taken to be a product of two simple factors Ni(c) and 
N 2 (o), and the basic iterative scheme can be expressed as; 

n n 

N (c) N (c) 6<t> = cuL<|) (12) 
1 2 

where o is a sequence of acceleration parameters, o is a 
relaxation parameter, Ni(a) and N 2 (a) are both functions of c, 
and they are easy to invert. The solution of Equation (12) is 
then obtained in two steps through an intermediate value F, that 
is 


n n 

step 1: N (a) F = ooL<}> 

1 

n n 

step 2: N(a)6<}> =F. 

2 


(13) 


This, in turn, requires an additional boundary conditions for F 
in order to solve for the step 1. Slow convergence or even 
divergence may occur if the values of acceleration parameters o 
and the boundary conditions for F are not carefully chosen. 
Although the effect of the intermediate boundary conditions in 
the AF scheme has been studied via the von Neumann and 
Gustaf sson-Kreiss-Sundstrom theory by South and Hafez^^, the 
stability analysis is valid only for purely subsonic flow 
calculations and there is still no rigorous analysis available 
for mixed subsonic-supersonic problems. Consequently, the 
performance of the AF scheme for transonic flow calculations may 
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be strongly depended upon the experience of an individual user. 
That is a fast convergence rate can be achieved if optimal values 
of acceleration parameters and suitable intermediate boundary 
conditions are provided, but not otherwise. 

Minimal Residual Algorithm 

In order to obtain a better approximation + i in the Newton 
like iterative scheme we need to solve a system of linear 
equations 

M 6<|) = - r (14) 
where M is a large sparse matrix operator. It is important to 
have an efficient solution method, since the linear system has to 
be solved for each step in the Newton like procedure. Direct 
method would not be possible since it requires a large amount of 
storages and arithemetic operations. In this paper an iterative 
method based on a Minimal Residual (MR) algorithm is used. 
Although the MR method has a slower convergence rate than the 
Conjugate Gradient algorithm^^, it can be applied to both 
symmetric and unsymmetric problems as long as all eigenvalues of 
the matrix operator has positive (or negative) real part. The 
number of iterations, NI , required to attain a given accuracy c 
using the MR method is given^^ by 

NI = 0.5 K ln(l/c) (15) 
where K= | |M| | | |M" ^ | | , is the condition number of the matrix 
operator M. Clearly the rate of convergence depends upon the 
value of K, in the sense that the smaller the value for K, the 
faster is the convergence rate that can be achieved. In order to 
accelerate the iterative process, a non-singular matrix C is 
introduced, and the linear system (14) is rewritten as 
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r 


(16) 


-1 

MC 6^ = - 


where 6^ = C64>. Equation (16) is known as the preconditioned 

system and C is the preconditioning operator. Suppose C is 

chosen so that C'^ is a good approximation to then the 

condition number of MC"^ would be much smaller than that for M 

itself. Consequently, solving the preconditioned system (16) 

will yield a faster convergence rate than that for the original 

system (14). A detailed account of the construction for C and 

its relationship to the matrix M will be given shortly. The 
% 

preconditioned MR algorithm can now be described as follows: 

Let 6<}>° be an initial guess correction vector, compute the 
residual vector, p“ = -r“ - M6(}>°, and solve Cz“ = p“, then for 
k=0 , 1 , 2 , . . . k , do 


k+1 k k 

6<|) = 6(f> + o z 

k 

k+1 k k 

p = p - a ,Mz 
k 

k+1 k+1 

sieve Cz = p 

k k 

where (p , Mz ) 


k k k 

(Mz , Mz ) 


(17) 


T 

Here (x,y) denotes the usual inner product, i.e. (x,y)=x y. The 
main computational work per iteration in the preconditioned MR 
algorithm is one matrix-vector multiplication for Mz, and the 
solution for Cz=p. 

Since we are interested in the overall convergence for the 


nonlinear problem (2), it may not necessary to solve the linear 
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system in (5) to excessively high accuracy for each Newton like 
iteration. In our implementation only an approximate solution is 
sought for each iterate, and this can be achieved by using a 
small fixed number of iterations (such as TT is set to 4) in the 
preconditioned MR algorithm. The iterative procedure described in 
this section thus consists of outer and inner iteration: the 
outer iterate is based on a Newton like algorithm (5), and the 
preconditioned MR method (17) is applied to find an approximate 
solution for the inner iteration. This method is therefore 
regarded as a Newton like - Minimal Residual algorithm, or simply 
an Inexact Newton like (IN) procedure. 

Preconditioning Matrix 

Having presented the preconditioned MR algorithm (17), we 
now study how to construct the preconditioning matrix C. 
Remember that the extra computational work for each iteration in 
the preconditioned algorithm is in solving the linear system 
Cz=p. Note that, if C=I, the identity matrix, then (17) becomes 
the regular MR algorithm. For a good preconditioned algorithm, C 
should be chosen so that the condition number, of MC‘^ is much 
smaller than that of M itself. In particular, if C=M the 
condition number of MC'^ is 1, the minimal value one could 
obtain. Although a converged solution can be achieved in one 
iteration, we need to solve Mz=p which is as difficult as solving 
the original problem M6((>=-r. Thus another important consideration 
is that C"^ should be easily invertable, otherwise the 
preconditioned algorithm will not be efficient. To satisfy these 
two criteria C is taken to be an approximation to the matrix 
operator M, and C is also a product of sparse triangular 
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matrices ; 


C = LU = M + E (18) 

where L and U are sparse lower and upper triangular matrices and 
E is the error matrix which measures how good is the 
approximation between C and M. The matrix C is based on an 
incomplete factorization technique^ ^ , and the algorithm for 
constructing the sparse matrices L and U can be described as 
follows. Recall that M is a sparse matrix consisting of five 
non-zero diagonals, where 

(M<t>) =E<}> +D(}) +F<}) (19) 

i/j i/j i,j i/j i-l,j i/j i+l,j 

+ H (f) + B <t> 

i/j i/j+1 i/j i/j-1 

Now L and U are constructed so that L and U has four non-zero 
diagonals respectively, in which the three non-zero diagonals are 
in the same positions to those in the lower and upper triangular 
part of M, where 

(Lcf)) = V (}) + t (}) 

i/j i/j i/j i/j i-l/j 

+ g <}) + x <|) (20) 

i/j i/j-1 i/j i+'l/j-l 

(U<j>) =<j> +e(}> +f<}> 

i/j i/j i/j i+l/j i/j i/j+1 

+ y <|> 

i,j i-l,j+l 

t 

The elements of L and U are computed from the coefficients of M 
according to the relations: 

g = B 

i/j i/j 

X = -g e 

i/j i/j i/j-1 
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t = D “ g y 

i/j i/j ifj i/j-1 

_ = (1 + a)E_ - g f 

i/j i,j i,j i,j-l 

-t _(e + y ) ( 21 ) 

i/D i-l.j i-l,j 

-X (y + e ) 

i,j 

e =(F -f X )/v 

i'li i/j i+l/j-1 i/j i,j 

y = -t f / V 

i/j i/j i-l/j i,j 

f = H /V 

i/j i/j i , j 

A small value , a, is added to the main diagonal elements of M 
to ensure the stability of the iterative scheme. The convergence 
rate, however, is not sensitive to o in the range from 0.025 to 
0.1, and o = 0.05 is used for all test problems in this paper. 
It should be pointed out that the algorithm given in (21) satisfy 
the following condition^®. 


non-zero elements 
(except the main 
diagonal) of M 


non-zero elements of C 
which are in the same 
locations as M 


Consequently although the preconditioning matrix C is factorized 

to LU decomposition, the product of LU will be symmetric as long 

as M is symmetric. Thus LU may be a symmetric matrix for purely 

subsonic flows or during the early iterations in mixed subsonic- 

supersonic calculations. However, when supersonic regions are 

developed a ^ term will apperar in the M matrix operator, and 
^t 

LU will then become asymmetric. 

The solution of Cp=z can now be obtained efficiently by the 
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following simple procedure. Since C=LU, the linear system Cz=p 
can be rewritten as 


Ls = 2 
Up = s 

where s is a dummy vector. The solution of Ls = 
through a forward substitution, where 


( 22 ) 

2 is obtained 


s - (z -g s -t s 

i,j i/j i/j i.j-1 i/3 i-1/3 

- X s ^ . 

i/j i+l/3"l i,j 

Note that, unlike the AF scheme, no boundary condition is 
required in order to solve for the dummy vector s. The solution 
of Up = s is obtained by a backward substitution, where 

p =s -e p -f P 

i,j i,j i,j i+l/j i,j i/j+1 

- y P 

i/j i-l,j+l 

To end this section we would like to mention that one of the 
differences in the implementation of the present method compared 
to that given in our early refenece"’ is that a more accurate 
preconditioning matrix is used in this paper. The computational 
results given in [7] were based on a different LU decomposition, 
in which x = y =0 for all i,j in Equations (20) (i.e. L 

i/j i/j 

and U had only three non-zero diagonals). With an additional 
diagonal for L and U respectively, the norm of the error matrix E 
is smaller than that in our early paper. Consequently, a faster 
convergence rate is achieved in the present iterative scheme. 
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IV. Numerical Results 

In this section results of numerical experiments using AF 
and IN iterative scheme are presented. The computer program is 
based on the TAIR code*, and they all have been carried out on 
the CDC CYBER 203 computer of NASA Langley Research Center. The 
problems to be considered are transonic potential flow fields 
around NACA 0012 airfoil at different Mach numbers and different 
angles of attack. The grid system used for both schemes is the 
same, where the total grid is 149*30 = 4470 points . Furthermore, 
the boundary conditions, initial starting vectors and criterion 
for convergence test are the same for both schemes. For all 
figures presented in this sections, the solid lines are the 
results obtained based on the AF scheme of Holst®, and the dotted 
lines are those for the IN iterative scheme. 

The surface pressure coefficient distributions, C , are 

P 

identical to those reported by Holst's experiments®, and hence 
they will not be presented here. In this paper, we shall mainly 
focus on the comparison of convergence rates and the efficiency 
for the AF and IN iterative schemes. Figures 1-3 compare the 
rates of convergence of the two methods for the following cases: 
(1) Moe = 0.85, c = 0"; (2) Moo = 0.8, a = 0.5“; and .(3) Moo = 0.75, 
a = 2®. The first one is for zero angle of attack, i.e. a 
non-lifting condition, (2) and (3) correspond to lifting airfoil 
calculations . 

From the comparison of convergence histories we observe that 
the IN iterative scheme produces generally a smoother reduction 
in the residual norm, especially for lifting airfoil 
calculations. Another useful criterion for comparing the 


17 


efficiency of each method is to study the development of 
supersonic points and the circulation as the number *of iterations 
is increased. Figures 4-5 show the development of the number of 
supersonic points for Moo=0,85, c=0" and Moe=0.8, o=0.5®, and 
Figure 6 gives the development of circulation for M»=0.8, 0=0.5". 
The results in Figures 4-6 clearly indicate that the number of 
supersonic points and the circulations are rapidly established 
for the IN iterative scheme. From Figures 4-5 we observe that at 
the 2nd iteration supersonic points had already been established 
in the IN scheme, and moreover, it reached almost 50% of its 
final value at the 5th iteration, whereas the AF scheme attained 
its first supersonic points in the 6th and 7th iterations for the 
cases of M«>=0.85, and Moo=0.8 respectively. 

It should be noted that the convergence histories alone do 
not reveal the complete picture for the comparison between the AF 
and IN schemes, another important point of concern is the 
computing time required for each method to obtain a given 
accuracy. For a detail comparison one should study the exact 
numbers of arithmetic operations for each computer code. 
However, this is obviously very difficult to achieve for 
practical problems, such as the transonic flow fields 
calculations. To provide a reasonable indication. Table 1 gives 
the total CPU time in seconds required to reduce an average 
residual (i.e. MrlU) to less than lO"’. The CPU time per 
iterate is 0.235 seconds for the AF scheme, and is 0.524 seconds 
for the IN scheme. 
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Table 1. Comparison of CPU 
time for the AF^na IN schemes 



Moo = 0.85 

GO 

O 

II 

8 

Mee = 0.75 


0 = 0“ 

c = 0.5“ 

0=2“ 

AF 

32.4 

35.5 

28.2 

IN 

41.8 

44.5 

31.3 


It should be pointed out that a considerable amount of 
computational work is needed for evaluating the residual vactor 
at each iterate, since it is necessary to update the fluid 
density at each grid point, modify the densities which are in the 
SjJpersonic regions, and calculate the velocity potentials, ... 
etc. In fact evaluating these residuals takes more computing 
time than solving the discrete potential equation using the AF 
scheme for each iterate. Consequently, although the work per 
iterate for the IN scheme takes twice the CPU time required by 
the AF scheme, the total number of iterations is reduced so that 
the overall computing time needed to attain the same accuracy for 
the IN scheme is not significantly larger than that based on the 
AF scheme. 

A further improvement in computing time for the IN iterative 
scheme is possible, since no effort to optimize the present 
computer code was attempted. One could expect that a larger 
number of inner iterations (i.e. the value for k in the 
preconditioned Minimal Residual algorithm) per outer interation 
(i.e. the value for n in the Newton like iterative process) might 
result in a smaller number of outer iterations, and similarly, a 
smaller number of inner iterations per outer iteration might 
result in a larger number of outer iterations. It is, of course, 
not clear what values are the best possible for a particular 
problem. However, a variable for the inner iteration, so that k 
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is gradually increased as the outer iterations proceed, will 
probably be a better choice than a constant value for k as used 
in the present implementation. A criterion for achieving this 
objective is being investigated. 

Finally, it should be mentioned that in assessing these 
numerical results, one should keep in mind that there is a major 
difference between these methods. Although the AF scheme 
requires less computing time, its formulation and application 
require specialized knowledge and experience of an individual 
user. In the sense that the performance of the AF scheme may be 
greatly affected by the choice of the acceleration parameters and 
also the treatment of the boundary conditions for the inter- 
mediate variable. Because of these reasons, many researchers 
have had a mixed experience with the AF scheme, sometimes finding 
excellent results, and sometimes finding them disappointing. 
Although more computing time is needed for the IN iterative 
scheme, it is easy to program and it does not suffer from the 
above difficulties, 

V. Conclusion 

A Newton like - Minimal Residual iterative scheme is pre- 
sented for the solution of the full potential' equation in 
transonic ranges. The method described here exhibited an 
attractive property over the Approximate Factorization scheme, 
namely the uncertainties and difficulties in choosing the 
acceleration parameters and treatment of boundary conditions for 
the intermediate variable are eliminated. Consequently, the 
present method is less problem dependent and also less user 
dependent as well. Numerical results for transonic airfoil 
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calculations are promising: the IN method produces generally a 

smoother reduction in the residual norm, * and the number of 
supersonic points and circulations are rapidly established as the 
number of iterations is increased. In addition, there is still 
room for improvement for the present method. Finally, two 
potential ares of application are as follows: 

(1) Transonic Wing Calculations: 

It is technically straightforward to extend the present 

method for numerical solution of a three-dimensional full 
potent ial • equat i on . Moreover, the increase of the computational 
work over a two-dimensional flow calculation is smaller for the 
present method compared to the corresponding increase for the 
Approximate Factorization scheme. For a three-dimensional 
problem, the matrix operator M in the Newton like iterative 
procedure will be a seven-point formula instead of a five-point 
as for a two-dimensional problem. However, the main computational 
work for the inner iteration is comparable to that required for a 
two-dimensional calculation, since a sparse LU factorization can 
be obtained with no difficulty. The Approximate Factorization 
scheme, on the other hand, consists of three-step calculations^*, 
and it can be expressed as 

n n 

N (o) N (c) N (a) 6<j> = ouLd (25) 

12 3 

The solution of this scheme would now require two additional 
boundary conditions for the two intermediate variables. 

(2) Finite Element Method in Fluid Dynamics Problems: 

Since the Approximate Factorization scheme is essentially 
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based on the alternating direction splitting methods, they will 
not be applicable since it is n<5 longer possible, in the finite 
element formations to partition the matrix operator in terms of 
the usual directional derivatives. The present method, however, 
does not suffer from this, restriction, since an incomplete sparse 
LU factorization can still be derived for the finite element 
matrix . 

These two areas of applications and other possibilities will 
be under investigation. 
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Figure Captions 


1. 

Comparison of 

convergence 

histories 

(NACA 0012 

airfoil , 


Moo=0.85, a=0°) 





2. 

Comparison of 

convergence 

histories 

(NACA 0012 

airfoil , 


M»=0.8, 0=0.5“ 

) 




3. 

Comparison of 

convergence 

histories 

(NACA 0012 

airfoil , 


Mt»=0.75, o=2“) 





4. 

Development of 

the number of 

supersonic 

points (NSP) as the 


number of iterations is 

increased 

(NACA 0012 

airfoil. 


Moo=0.85, 0 = 0 °) 





5. 

Development of 

the number of 

supersonic 

points (NSP) as the 


number of iterations is increased (NACA 0012' airfoil, M<»=0.8, 
c=0.5“ ) 

6. Development of the circulation (CL) as the number of 
iterations is increased (NACA 0012 airfoil, Moo=0.8, a=0,5“) 
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